Exact phase diagram of quasispecies model with mutation rate modifier 
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We consider an infinite asexual population with a mutator allele which can elevate mutation rates. 
With probability /, a transition from nonmutator to mutator state occurs but the reverse transition 
is forbidden. We find that at / = 0, the population is in the state with minimum mutation rate and 
at / = f c , a phase transition occurs between a mixed phase with both nonmutators and mutators 
and a pure mutator phase. We calculate the critical probability f c and the total mutator fraction 
Q in the mixed phase exactly. Our predictions for Q are in agreement with those seen in microbial 
populations in static environments. 

PACS numbers: 87.10.-e, 87.10. Ca 



In the absence of recombination, biological evolution is 
driven by two competing processes namely selection that 
tends to localise the population around a fitness peak 
and mutation which has the opposite effect of delocalising 
it. Extensive theoretical and experimental studies have 
shown that there exists an error threshold beyond which 
the mutational load (fitness reduction) becomes too high 
to be compensated by selection pressure P, Q • For this 
reason, and because most mutations are known to have 
deleterious effect 0, 0, H, 0] > the spontaneous mutation 
rate is expected to be minimum subject to physicochem- 
ical constraints and physiological costs 0, 0] • 

However, mutators with mutation rate higher than the 
wild type can be produced when for example, an organ- 
ism is unable to neutralize mutagens such as radiation 
or repair DNA damage during replication In fact, 
hypermutable strains in high frequency have been found 
in some cancerous cells [10( and in natural isolates of 
certain pathogenic bacteria which persist for many years 
inspite of the presence of antibiotics Subpopula- 
tions with 10— 100 times higher mutation rates have also 
been seen to arise spontaneousl y in long term adapta- 
tion experiments on E. Coli [II 111 EI 111, El]. Due to 
the persistence of mutators at long times, it is important 
to study their role when mutation-selection balance has 
been reached. 

Here we consider the mutator problem within the 
framework of quasispecies model which is defined on the 
genotypic space and assumes an infinite population 
The transition from a nonmutator to mutator state oc- 
curs with probability / but the reverse reaction is ignored 
as it has a much smaller probability than / 5]. While 
the mutators are thus continually generated, they are 
selected against due to high mutational load. At suffi- 
ciently high /, we may expect the mutators to take over 
the whole population. However the possibility of such 
a transition has not been considered in previous stud- 
ies on smooth fitness landscapes 17, El, 111, HI]- Here 
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FIG. 1: Schematic phase diagram of the quasispecies model 
with mutation rate modifier. With probability /, the nonmu- 
tators (mutation probability /i) change to mutators (muta- 
tion probability v) . The pure nonmutator phase occurs when 
/ = and pure mutator phase for / > f c . The system is in 
the mixed phase for < / < f c . 



we show that the system can be in one of the following 
three phases (see Fig. Q]): a phase with only nonmuta- 
tors (/ = 0), a mixed phase in which both mutators and 
nonmutators are present (/ < f c ) and a phase with only 
mutators (/ > f c ). The critical probability f c at which 
a phase transition occurs between the mixed phase and 
a pure mutator phase is found exactly. We also calculate 
the total mutator fraction exactly for which approximate 
expressions have been obtained in the past [19|, l20| . 

We consider a haploid asexual population of infinite 
size evolving in a static environment. An individual 
in the population is represented by a binary sequence 
a = {a 1 , er^} with L loci where a. L — or 1. Each 
sequence is endowed with fitness F(a) which is propor- 
tional to the average number of offspring produced per 
generation. As there is considerable experimental evi- 
dence that the individual loci can contribute indepen- 
dently to the genome fitness 2l|, we work with multi- 
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plicative fitness F(a) = - =1 (1 — s)"" 1 where the selec- 
tion coefficient s G [0,1). Assuming that independent 
point mutations occur during the replication process, a 
sequence a' mutates to a with a probability given by 
the mutation matrix M M (<7, a') = ^ d O' ff ')(i - ^) L - d (^^') 
where d(<r, a') — $^i=i &i + o\ — 2(7^ is the Hamming 
distance between the sequences a and a' and fi is the 
mutation probability per locus per generation. The mu- 
tation rate modifier allele is modeled by attaching an 
additional bit with each sequence which controls the mu- 
tation rate but does not affect the fitness. A nonmutator 
sequence has a spontaneous mutation probability /i per 
locus per generation while the one with mutator allele 
mutates with probability v = A/z where A > 1. With for- 
ward probability /, the mutator allele is obtained from 
the nonmutator and the reverse reaction is neglected. 

Thus the average fraction P(a,t) and Q{a,t) of the 
nonmutator and the mutator respectively at generation t 
evolves deterministically according to the following cou- 
pled nonlinear difference equations: 
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The average fitness W(t) = £ CT F(V) [P(a, t) + Q(a, t)] 
in the denominator of the above equations ensures the 
normalisation condition P(a, t) + Q(a,t) — 1 is sat- 
isfied. We are interested in the steady state when these 
frequencies become time-independent. The steady state 
phase diagram is summarised in Fig. [TJ If / = 0, al- 
though the steady state fractions P{cr) and Q(cr) obey 
similar equations, the minimum mutation rate is chosen 
and the population is in a pure nonmutator phase with 
a quasispecies localised around the fitness peak as there 
is no error threshold for multiplicative fitness [H- For 
/ 7^ 0, while the nonmutator population reduces due to 
nonzero forward rate, they are favored over the muta- 
tors since the latter have the tendency to delocalise due 
to elevated mutation rates. Due to this competition, we 
may anticipate a phase transition in the A — / plane be- 
tween the mixed and pure mutator phase. Note that this 
transition is different from the error threshold transition 
in which the population delocalises from the fitness peak 
beyond a critical error rate 0. 

Before discussing the phase transition, we first demon- 
strate that for generic initial conditions, the popu- 
lation is in the state with minimum mutation rate 
when / = 0. Writing P(a,t) = Y(a,t)/X(t) and 
Q(a,t) = Z(a,t)/X(t) in Eqs. Q and © where X(t) = 
Y, a ,Y{p' ) t)+Z[p',t) and 



we find that the unnormalised variables Y(a, t) and 
Z(a, t) obey linear uncoupled equations given by [2| 

Y(a, t + l)=J2 M ^ <?')F(a')Y(a', t) (4) 
Z(a, t + l) = J2 M "( CT > v')F{a-')Z{o-' , t) . (5) 

The solution of the above equations is of the following 
product form, 

L L 

Y(a,t) = Ylvl^^y? > z (v,t) = n z l~ a * z ? (6) 

i=i i=i 

where the time-dependent fractions yk and Zk , k = 0, 1 
obey Eqs. ^ and for L = 1. This can be verified, for 
example, for Y(a, t) by using the above ansatz in Eq. (|4|) 
whose right hand side (RHS) can be expressed as a prod- 
uct over L terms, 



a' i=l 
L 



1-fi 



= II ^ 0- - rf~ ai yo(t) + M 1 ^' (1 - ^ (1 - s)yi (t) 

i=l 
L 

= IIfo""*(* + i) yi € (* + i) 



where we have used the evolution equations for yo,yi 
to arrive at the desired result. For an initial condition in 
which all the population is at the least fit sequence — 
{1,1, ...,1} with unnormalised nonmutator population 
a 7^ and mutator population j3, the one locus frac- 
tions yk and Zk can be straightforwardly computed and 

( \ 1 / L t i \ t i \ 

where 



x(t + 1) = F (°) t) + z & *)] 



(3) 



(2- S )(l- M )± ^(1-^)2 
K± = 2 ' (?) 

Since k_ < k + and n + (y) < it follows that z /y 

vanishes when t — > oo. Using this in the expression for 
the average fitness W(t), it follows that the steady state 
fitness W = K+([i) as in the pure nonmutator phase. 

We now consider the interesting situation when the 
forward rate / is nonzero. In the following, we will first 
find the average fitness W < for f < f c and W> for / > f c 
and then determine the critical point f c by matching W< 
and W> at the transition. The steady state equations 
for the population fractions, unlike the time-dependent 
ones, can not be linearised by passing to Y and Z vari- 
ables due to Eq. ((3]). As a consequence, due to the nor- 
malisation factor W in the denominator, Eq. (p} for the 
nonmutator fraction is coupled to the mutator fraction. 
However on summing over a on both sides of Eq. |1} in 
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FIG. 2: Average fitness W as a function of / (main) and A 
(inset) for L = 20, fi = 10~ 4 and s = 0.05. 



the steady state, we find that provided the nonmutator 
fraction P = P(a) is nonzero, the average fitness W 
does not depend on the mutator fraction and writeable 

as 



W = 



(i-/)£ g W» 



P^Q 



(8) 



thus leading to an uncoupled nonlinear equation for P(<r). 
Eliminating W from Eq. (JTJ) using the above equation and 
writing P(a) — P (c)/£ ff ' P {<?'), we see that P{a) obeys 
the quasispecies equation for a population without mu- 
tation rate modifier. In this case, the normalised steady 
state distribution is known to be given by [22j , 



(9) 



where f>o,pi are the solutions of the corresponding one 
locus model and the average fitness is given by 
as seen in the pure nonmutator phase. From Eq. ([5]). we 
thus obtain W< = (I- f)J2 a F(a)P(a) = (1-/) K$(ji). 
Contrary to naive expectation, the fitness W < is unaf- 
fected by the mutation rate v. This result is consistent 
with the reduction principle that requires the average fit- 
ness to be maximised [7], [8[ . This can happen if only the 
nonmutators contribute to the fitness thus minimising 
the mutational load due to mutators. But as the forward 
rate is nonzero, this contribution is reduced by a factor 
1 — /. Besides P ^ (mixed phase), P — is also a so- 
lution of Eq. . For f > f c phase in which this solution 
is valid, the mutator population Q(a) = Yif^i^cT^'^? 



(see Fig. [3]) and the average fitness W> — 

A plot of average fitness as a function of / and A is 
shown in Fig. With increasing / as the population 
goes from the mixed to pure mutator phase, the average 
fitness decreases to a constant since the nonmutator frac- 
tion vanishes for / > f c . As shown in the inset, in the 
pure mutator phase, the fitness decreases with increasing 
A since the increase in mutation rate decreases the pop- 
ulation at the top of the fitness landscape. In the mixed 



FIG. 3: Phase diagram in the / — A plane obtained using 
Eq. (TTDJ with fiL = 0.003 and s = 0.05. The horizontal lines 
show the maximum possible value of A as v = A/i < 1. Inset: 
Mutator distribution Q(d) as a function of the distance d from 
the master sequence for L = 10 and f c — 0.025. 



phase, as discussed above, the fitness is constant in A. 
Matching the fitnesses W < and W > at the critical point, 
the phase boundary in the / — A plane is obtained (see 
Fig.©, 



(l-/c) 



1/L 



(2 - s)(l - U C ) + y/4^(l - s) + s 2 (l - v c f 

(2 - - n) + vva-^ + ^a-A*) 2 

(10) 

When v = or s = 0, the critical point f c = and the 
population is always in the pure mutator phase. 

So far we have discussed the quasispecies model for 
arbitrary genome length L. To compute the total frac- 
tion Q = Q(a) of mutators in the mixed phase, we 
now consider the model defined by Eqs. ([T]) and ([2]) when 
the genome length L — > oo and the mutation probabil- 
ities fx, v —* with U — fiL and V — vL finite. In 
this limit, a sequence at a Hamming distance k from the 
fittest mutates to one at distance j with mutation proba- 
bility Mu(J, k) which is a Poisson distribution with mean 
U for k < j and zero otherwise [22L I23I ] . Furthermore, 
the average fitness W < = (1 — f )e~ u and M / > = e~ v so 
that the critical probability f c = l — e~ v+u , independent 
of s. The fractions P(k) and Q(k) of nonmutators and 
mutators respectively at a Hamming distance k from the 
fittest sequence obey the following equations [20| . 



P(k) 



Q(k) 
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(11) 
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fe'=0 v 



k-k' 



k') 



il-sfQ(k') + c 2 P{k)(12) 



where the coefficient c\ = (1 — / c )/(l — /) and ci — 
f /(I — f). Using the Eqs. (fTT|) and (fT2|) . one can write 
down the recursion relations for the generating function 

G(z) = EZo z k Q(V and H(z) = ZZo zkp ( k ) wi * h 
G(l) = 1 — H(l) = Q. On applying these recursion 
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FIG. 4: Variation of mutator fraction Q with s (main) and 
A (inset) in the mixed phase obtained using Eq. (If 3[) when 
n — » oo. 



equations n times and setting z — 1, we obtain 

Q = c5 l e y ^=« (1 ~ s) ' i G'((I-s)™) 



+ ca(l-Q) 



(13) 



m=0 



In the limit n — > oo, the first term on the RHS drops 
out for ci < 1 (/ < fc) and the sum can be expressed 
in terms of an incomplete gamma function by replacing 
the infinite sum by an integral. For biologically realistic 
situations for which / «s«F-(/ (see below), the sum 
can be calculated approximately to yield 



Q 



fV27T 



fV2t+(i-fW(y-u)s 



(14) 



which increases with / but decreases with U, A and s (also 
see Fig. 0}. 

As an application of our results, we consider a large 
population of bacteria E. Coli for which the genome mu- 
tation rate U « 2 x 10 -4 0, H, [l6| and the spontaneous 
forward transition rate / w 5 X 10~ 7 [24| has been es- 
timated. Our calculations show that a subpopulation of 
weak mutators with A = 10 can take over the entire popu- 
lation if / > 2 x 10~ 3 . Such a high transition rate can for 
example occur in the presence of mutagens (l4j |. More- 
over, our preliminary simulations on finite populations 
indicate that f c is considerably reduced due to stochas- 
tic fluctuations and it should be possible to observe this 
phase transition in experiments even without mutagens. 
As shown in Fig. |4j we obtain Q ~ 0.5% - 0.03% for 
s ~ 10~ 5 - lO" 1 and A = 10 using Eq. ([13]). Such small 
fractions have been observed in several long term exper- 
iments on E. Coli: 0.25% in 



12j, 0.6% in [13j and more 



recently, 0.5% in [14|. 

To summarise, we have presented several exact results 
for a quasispecies model with mutation rate modifier. 



The model discussed above can well describe large pop- 
ulations in a stable environment and harboring less than 
1% of mutator subpopulation . However if the popula- 
tion is exposed to a continuously changing environment 
as in infectious diseases 11] , the mutator fraction can be 
as high as 50 — 70% and we need to consider the evolution 
on dynamic fitness landscapes [IH, 03 • 
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